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ABSTRACT 


This thesis formulates the full nonlinear equations of motion for determining the 
stability of helicopter coupled rotor-fuselage motion utilizing MATLAB®‘s Symbolic 
Math Toolbox. Using the extended symbolic processor toolbox, the goal of this work was 
to eliminate the time consuming process of converting Fortran or C code generated by the 
symbolic processor, MAPLE® into a MATLAB® useable format where it is further 
incorporated into an ‘S-function’ to be used in the dynamic simulation environment. 

The formulation of the equations of motion utilized in this process is unique in 
that it uses the complete set of nonlinear terms in the equations of motions without 
utilizing ordering schemes, small angle assumptions, linearizing techniques, or other 
simplifying assumptions. After derivation, the equations of motion are numerically 
integrated using the dynamic simulation software SIMULINK® and a time history plot is 
generated of blade and fuselage motion. The equations of motion are regenerated with 
each time step allowing the adjustment of characteristic structural, blade and dampening 
properties. These time traces can be used to explore the effects of damping 
nonlinearities, structural nonlinearities, active control, individual blade control, and 


damper failure on ground resonance. 
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I. INTRODUCTION 


A. DISCUSSION 

Air and ground resonance are potentially destructive mechanical instabilities that 
may occur in helicopters due to dynamically coupled interaction between the rigid body 
motion of the fuselage and oscillations of the rotor blades in their plane of rotation. For 
over 40 years engineers have felt reasonably comfortable with their ability to predict and 
control air and ground mechanical instabilities. The self excited vibrations that occur 
while an aircraft is : contact with the ground is known as ground resonance and this will 
be the primary focus of this thesis. 

Helicopter manufacturers employ several methods of blade root retention to 
obiainthe flapping and lead-lag degrees of freedom necessary for a helicopter rotor blade 
to operate properly. These methods include: (1) fully articulated; (2) bearingless; (3) and 
hingeless main rotor designs. 

Fully articulated rotor heads utilize two hinges that rotate with and are attached to 
each rotor blade in order to provide the necessary freedom of movement for the dynamic 
motion of the rotor blade. One hinge is positioned in a horizontal plane and allows for 
flapping motion of the rotor blade. The other hinge is positioned in a vertical plane. This 
hinge allows for lead-lag, movement of the blade in its plane of rotation. The lead-lag 
hinge generally incorporates some damping. This is generally a hydraulic air-oil strut or 
an elastomeric design, to constrain the movement of the individual rotor blades. Early 


designs of these dampers have tended to be heavy and maintenance intensive. 





The hingeless/bearingless rotor is a variation of the fully articulated rotor system. 
In this configuration the designer has substituted a flexible section to replace both the 
lead-lag hinge as well as the flapping hinge. 

In short, by eliminating the need for hinges, improvements in technology, 
particularly composite materials, have allowed designers to significantly reduce the 
weight and complexity of modern helicopter rotor systems. However, this 
eiccon/eliiination of hinges has also reduced the likely attachment points for linear 
dampers and other means of controlling the motion of the rotor blades. This further 
complicates eliminating ground resonance, as well as the associated complexity of 
manufacturing and maintenance of the dampers. By simplifying and facilitating the 
process of modeling ground resonance the hope is to open the study of ground resonance 
to a deeper understanding as well as investigate methods of controlling it. By allowing 
the variation of parameters and by receiving immediate feedback as to the resultant 
effects, this program and future versions will help improve the properties of existing 
damper systems, facilitate the possibility for application of nonlinear dampers, potentially 
eliminate external dampers using existing material technology, and facilitate the 


investigation of other aero-mechanical instabilities. 


B. BACKGROUND 


In the early 1940’s the phenomenon of ground resonance was already being 
investigated by NACA. Robert Coleman’s pioneering work was the first definitive 
analysis to correctly address the issue of ground resonance [Ref. 1] and [Ref. 2]. In his 


work Coleman identified ground resonance to be a self-excited, mechanical instability 


phenomenon. He found the primary modes to be excited in the normal operation of the 
rotorcraft were the hinged rigid body response of the blades, in their plane of rotation 
coupled with horizontal deflection of the main rotor pylon. 

M. L. Deutsch simplified the results of Coleman and Feingold for the practicing 
engineer. Based on Coleman’s theory he developed a quantitative analysis of the 
damping required to keep the helicopter free of the instability, often referred to as 
Deutsch’s Criteria [Ref. 3]. 

Coleman and Feingold’s work became the basis for the evolution of theory and 
design techniques used for dealing with ground resonance. Although the Siesac theory 
offers much insight and understanding into the phenomenon, especially for conventional 
articulated rotor systems, the increasing popularity of hingeless and bearingless rotor 
designs in modern helicopters called for increasingly more sophisticated analytical tools. 

More sophisticated analytical tools became possible with the evolution of digital 
computers and the improvement in computational power. Peters and Hohenemser [Ref. 
24] apply Floquet analysis to the problem of lifting rotor stability. Floquet analysis is a 
method which can be used to determine the stability of solutions to systems of linear 
ordinary differential equations with periodic coefficients. The Floquet transition matrix 
which relates the system state variables at the beginning and end of a rotational period is 
computed by numerical time wise integration. The eigehvalund of the transition matrix 
are a measure of system stability. Hammond [Ref. 25] applies Floquet analysis to the 
prediction of mechanical instabilities, specifically examining the case of unbalanced lead- 


lag damping. The unbalanced problem requires solution of the equations of motion with 





the blade equations expressed in the rotating reference frame because a transformation to 
the fixed system is no longer possible for a ground resonance analysis as was possible for 
the isometric case. As a result, you are left with a system of equations with periodic 
coefficients which can be handled by the Floquet method. 

Hingeless and bearingless rotor configurations often face the additional difficulty 
of air resonance. Aerodynamics may play more than a passive roll in the ground 
resonance regime in hingeless systems in contrast to articulated systems where 
aerodynamics have little effect. As a result, more complex models are required to 
accurately represent the physics of the helicopter aeromechanical stability problem. 
Models must include blade flap and torsional degrees of freedom as well as lead-lag 
degrees of freedom. Fuselage models also should include pitch and roll as well as 
translational degrees of freedom. Aerodynamic models can range from quasi-steady strip 
theory to unsteady aerodynamic theories which include elaborate wake models or 
dynamic inflow models. Ormiston [Ref. 26] utilizes a rigid blade and rigid fuselage 
model with flap-lag and pitch-roll pikes of freedom to conduct parametric 
investigations based on an eigenvalue analysis. As is typical, the equations of motion 
were derived by a Newtonian approach and the resulting system of nonlinear differential 
equations are linearized for small perturbations. The model includes linear rotor blade 
and landing gear springs, viscous damping, and quasi-steady aerodynamics. Freidmann 
and Venkatesan [Ref. 27] and Freidmann and Warmbrodt [Ref. 14] derive the complete 
set of governing equations of a helicopter rotor coupled to a rigid body fuselage. The 


equations account for rotor blade elastic deformations and include quasi-steady 





aerodynamics or modified Theordorsen unsteady aerodynamic theory. In deriving the 
full equations of motion, Freidmann et al., stress the importance of applying an ordering 
scheme to the process in order to handle the complexity of the equations and enormous 
number of terms generated by their expansion. The equations, as presented by Freidmann 
et al. [Ref. 27 and 14], are in a form which makes them generally applicable to a wide 
range of rotorcraft problems. | 

Another interest in the study of helicopter ground resonance is the effect that 
nonlinear elastic and damping forces have on stability. Tongue [Ref. 23], Tongue and 
Flowers [Ref. 9 and 8], Tongue and Jankowski [Ref. 28], and Tang and Dowell [Ref. 29], 
use variations of the nonlinear technique of harmonic balance using describing functions 
to represent nonlinear damping. The technique is useful for investigating limit cycle 
behavior of strongly nonlinear systems and its impact on system stability. 

Helicopter aeromechanical instabilities can be analyzed by methods ranging from 
Coleman’s classic analysis to direct time integration of the equations of motion. As . 
engineers strive to develop rotor systems free of ground and air resonance which do not 
require the addition of maintenance intensive mechanical damping systems, more 
elaborate models will be needed to accurately capture all physical aspects of the problem. 
To achieve the truly damperless rotor Ormiston [Ref. 30] addresses three different 
approaches which may be feasible, 1) incorporating high damping material into the blade 
or flexbeam structure, 2) automatic feedback control, and 3) Pen ere of aeroelastic 
couplings to provide inherent stability. These three approaches provided the impetus 


behind the work performed by LT Christopher S. Robinson. 








A full nonlinear simulation model of coupled rotor/fuselage interaction was 
created by Robinson in March of 1997 at the Naval Postgraduate School (NPS) [Ref. 4]. 
This model is a dedicated Coleman analysis tool that initially utilizes a symbolic 
processor, Maple® to derive the full nonlinear equations of motion of a helicopter using 
the LaGrange equation. In Robinson’s thesis the derived equations of motion are then 
converted, in Maple®, to Fortran or C and then converted into a MATLAB® 
programming language. The converted MATLAB® result is incorporated into a 
SIMULINK® S-function, producing time history plots of blade/fuselage motion. This 
process is unique in that it uses the complete set of nonlinear terms in the equations of 
motions without utilizing ordering schemes, small angle assumptions, linearizing 
techniques, or other simplifying assumptions. 

Further development of the NPS modeler was accomplished by Robinson with the 
help of LCDR Robert L. King, [Ref. 12, 18, 19, 20]. The modeler was used to simulate 
the Froude Scale RAH-66 [Ref. 21]. It was also used to investigate an interblade 
coupling approach to improved rotor stability without lag dampers [Ref. 22]. 

King completed his dissertation in June 1999 [Ref. 5]. In his dissertation, King 
adopted Robinson’s work to further explore the potential of eliminating the snubber- 
aatiped or damper on hingeless rotor designs and replace it with a flexbeam that has been 
modified to possess nonlinear properties. It is King’s research in this field that has 
provided the motivation to accomplish the following work. 

Since Robinson’s original work, MATLAB?® has incorporated a symbolic 


processing toolbox into its software. By using the extended symbolic processor toolbox 





in MATLAB® the goal is to eliminate the awkward and time consuming process of 
converting the Fortran or C code generated by MAPLE® into a MATLAB® useable 
format where it is further incorporated into an ‘S-function’ to be used in the dynamic 
simulation environment. Direct interface between the derivation of the equations of 
motion and the actual simulation process results in fewer steps required to develop a 
system, it reduces the computer programs required to run a simulation, and reduces the 
potential inaminsesh of the solution that may result from human interaction. 
Hopefully this will encourage further use of the rotor-fuselage coupled motion simulation 
and enable researchers of all levels to explore the phenomenon of ground resonance. 
After a brief explanation of the work performed by Robinson, Rafanello, and 
King an introduction to MATLAB®, MATLAB® the symbolic processing toolbox, the S- 
function, and the Simulink model will be provided. The goal, again, is to provide a path 


for future use of the NPS modeler. 
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Il. CURRENT MODEL 


A. REVIEW OF WORK PERFORMED BY ROBINSON 


1. Simple Mode] 

As with this thesis, Robinson’s work began with a simplified model using reduced 
degrees of freedom on a three bladed model. The simplified model is based on that of 
Coleman as is shown in Figure 1, [Ref. 1]. Elastic forces generated by rotor blade and 
flexbeam motion were modeled as a linear torsional spring located at the effective hinge 
position of the blade. Landing gear stiffness was also modeled with linear springs. The 
landing gear and lead-lag dampers were modeled with linear dashpot type dampers. 
Transformation between the various systems, in order to develop Lagrange’s equations of 
motion, were accomplished using Euler angle rotations. For the simple model the 
coordinate transformations used were hub (parallel to fuselage but offset a distance h in 
the positive z direction) to inertial (fixed relative to the earth), blade undeformed (fixed to 
the effective hinge position) to hub, and blade deformed (fixed to the effective hinge 


position with the x-axis coincident with the blade ‘elastic’ axis) to blade undeformed. 








Damper 


5 AR, AA 








Figure 1. Simplified Rotor Model [From Ref. 4] 


Robinson showed excellent agreement between his model and the Coleman 
model, with departure occurring only when displacements are very large, Figure 2. This 
is to be expected since Robinson’s model does not assume smal] angle theory whereas the 


Coleman-Feingold-Bramwell model does. 
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Comparison of Simulation Results to Solution of Coleman Equations 
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Figure 2. Validation of model with solution of Coleman’s model [From Ref. 4] 


Figure 3 is a parametric plot of the simple model with an isotropic pylon and 
rotor. The damping ratio from the moving block result is plotted versus 0)/«, for various 


Deutsch numbers [Ref. 15]. 
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Moving Block Results Parametized with Deutsch Criteria 
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Figure 3. Validation of model against Deutsch Criteria [After Ref. 4] 


This model allows for the following degrees of freedom: 


ul - Fuselage translation in 1-direction (x-direction) 


u2 - Fuselage translation in the 2-direction (y-direction) 


Cx- Lead-lag angular displacement of k™ rotor blade 


2. 


Complex Model 


The complex model is based on that used by Straub [Ref. 6]. This model assumes 


rigid blades and fuselage, as in the simple model, with flap and lead-lag hinges as 


coincident. It incorporates fuselage rotation as well as blade flap, requiring additional 


12 


transformations. Figures 4 and 5 show a schematic representation of the complex model 
and the transformations utilized in the derivation. The transformations required in the 
complex model are fuselage (fixed to center of gravity of fuselage) to inertial, hub to 


fuselage, blade undeformed to hub, blade deformed to blade undeformed. 
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Figure 4. Schematic representation of complex model transformations 
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Figure 5. Schematic representation of complex model transformations 


This model utilizes the following degrees of freedom: 
ul - Fuselage translation in 1-direction (x-direction) 
u2 - Fuselage translation in the 2-direction (y-direction) 
rl — Fuselage rotation about 1-axis (roll) 
2 ~— Fuselage mere about 2-axis (pitch) 
Cx- Lead-lag angular displacement of k" rotor blade 


‘Bx - Flap angular displacement of k" rotor blade 
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3. Derivation of the Lagrange Equation 


A brief overview of the derivation accomplished by Robinson will be included in 
this thesis. For a more extensive explanation see Robinson’s thesis [Ref. 4]. The 


Lagrangian equation may be expressed as follows: 


Where T is the kinetic energy, U is the potential energy, D is the dissipation 
function, F; is the generalized force, and q; is the generalized displacement. In the 
complex model the generalized force term Fj, will describe the aerodynamic forces on the 
individual rotor blades. This derivation develops a system of homogeneous equations. 
The various energy terms are broken down into two categories, blade motion and 


fuselage motion to give the following equations: 


N 
T=T, +>), 
k=l 
N 
U=U,;+>, Us) 
k=] 


D=D,+>\(Us5); 


k=] 


Where the subscripts F and B indicate the fuselage and rotor blade respectively. 


a. Kinetic Energy Terms 


The kinetic energy of the kth rotor blade is given by the following 


expression: 
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ay a & 
(Tq) = |=m'@B- par 
0 


Here / is the position of a point on the elastic axis of the k™ rotor blade 
with respect to the inertial coordinate system at any instant in time, and m’ is the mass 
distribution per unit length of the blade (mass distribution per unit length is assumed to 
be uniform). The position of a point on the elastic axis of a rotor blade, 0, is expressed 
as the sum of relative positions with respect to the various coordinate systems 


transformed to the inertial system. Thus, 
p = (Or 1 ); a5 (Pz ¢ )r BG (Osy_ x )r + (Peg py )r + (Pp_pp)s ? 
where, for example, the term ((,, ;,), is the position of the origin of the 


undeformed blade coordinate system with respect to the hub coordinate system 
transformed into the inertial coordinate system. 

. Following the proper transformations, the resultant expression provides 
the position of an arbitrary point on the elastic axis of the k" rotor blade, with respect to 
the inertial coordinate system, at any instant in time, in terms of the system degrees of 
freedom. The time derivative of this expression gives the velocity, p , which is 
substituted into the equation for kinetic energy for the k™ rotor blade. 

The fuselage kinetic energy in terms of translational and rotational degrees 
of freedom is 


1 <i. hl ; 
(T; joo = 5M ait + Malls 


Dips esorss fhe be a 
(Te yor = raul +5 lah — 21 hh 


b. Potential Energy Terms 


The elastic forces generated by rotor blade motion give rise to a potential 
energy term in the Lagrange equation. Since a rigid blade model was assumed, the 


potential energy was modeled using equivalent torsional springs to restrain the rotor 
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blade, with spring constants selected to approximate elastic forces due to in plane and out 
of plane bending of the rotor blade (and the flexbeam in the hingeless case). The 


potential energy of the kth rotor blade is 
sul a eee 
(U;), =3 Koh ra 


The fuselage potential energy in terms of translational and rotation degrees 
of freedom is 
ene eee | 
(UU; 1 ee =e +5 Kas 
1 21 2 
UF) ror => KR + —KR,r, 
2 2 
An explanation of the validity of using an equivalent torsional spring 
system to model the elastic forces of a deformed rotor blade is given in some detail in 


Venkatesan and Friedmann [Ref. 7]. 


c. Dissipation Function Terms 


System damping is modeled in energy form by use of a dissipation 
function, which for the kth rotor blade of he complex rotor model is 
(Dy), = 4,82 +2.0,82 
2 2 
The dissipation function for the fuselage in terms of translational and rotational 
degrees of freedom is 
1 1 


(D; Jig 2 Cu; + 2 Cy; 


(De ) ror =>C Ri +5C Ry, 


d. Resultant Equations 


It is important to note that when applying Lagrange’s equation in 


MAPLE®, derivatives with respect to the degrees of freedom and the time rates of change 
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of the degrees of freedom, the time functional notation which represents these variables 


must be converted to independent variable notation. For example, the flap angle degree 


of(t) 
ot 





of freedom, B(t), and its time rate of change, , would have to be replaced in all of 


_ the energy expressions by the independent variables, B and dB respectively, in order for 


3a; er 


terms such as 0g; and (where q; = B(t) and or ) to be evaluated properly by 


alas, 

the MAPLE® symbolic engine. Additionally, for the time derivative term, or\ 04, , to 
be evaluated properly, all degrees of freedom expressed in independent notation must be 
converted back to tine dependent notation. 

The equations of motion generated for the simplified model (Coleman 
model) in this format were compared, by Robinson, to the equations used by Flowers and 
Tongue [Ref. 8,9]. The equations were found to match exactly except for the nonlinear 
damping terms that were not included initially but are later included in the simulation. 
Robinson’s work has been further utilized by LT Salvator Rafanello and LCDR Robert L. 


King. 


B. REVIEW OF WORK PERFORMED BY RAFANELLO 


Rafanello’s work further validated the NPS ground resonance modeler by 
matching simulation results with Professor Roberts E. Wood’s HSS-2 ground resonance 
analysis [Ref. 10, 11]. Rafanello used a five bladed version of the simplified model 


adapted to the H-3 Sea King. A mass-spring-damper system of the H-3’s landing gear 
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was examined to obtain parameters to enter into the NPS simulation. A thorough 
analysis of the natural frequencies in the lateral and roll modes at three power settings of 
0/20/80 percent airborne were made with respect to the literature prepared by Coleman, 
Feingold, and Deutsch. Stability charts were developed from the tailored data of the Sea 
King and compared with the modeler’s results. Through these results Rafanello was able 
to directly link the classical work with the NPS modeler in the roll and lateral modes of 
the H-3. 

In his work, Rafanello was also able to verify the conservative nature of 
Deutsch’s Criteria as was presented in Robinson’s thesis. He was able to determine that | 
Deutsch’s Criteria is conservative in that the critical point for his criteria, or bucket of 
each curve, still shows positive damping when analyzed with the NPS simulation. Please 


refer to Reference 10 for further analysis. 


C. REVIEW OF WORK PERFORMED BY KING 


King’s work focused on examining alternative designs for a damperless helicopter 
rotor. He explored the potential of eliminating the snubber-damper or damper on 
hingeless rotor designs and replacing it with a flexbeam that is modified to possess 
nonlinear properties. The NPS ground resonance modeler is well suited to this task in 
that it can accurately model nonlinear mechanical properties. King included nonlinear 
properties in the blade root to produce potentially acceptable bounded responses in the 
parameter region, with the NPS modeler, where linear theory would have predicted 
instability. He was also able to determine regions of unacceptable response where linear 


theory would have predicted stability. 
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King also looked at linearly linked rotor systems and unevenly spaced rotors 





using the NPS modeler. By linearly linking the rotors he was able to change the mode 
shapes of the rotor blades’ lead/lag motion. This alters the blade natural frequency in the 
lead/lag and regressing mode of the rotor. Therefore, by linking the blades, the potential 
exists to avoid ground/air resonance by detuning the rotor-body dynamics. This method 
avoids ground and air resonance in a similar manner to the nonlinear stiffness approach. 
Unevenly spacing the rotor blades produced the same mechanical instabilities as were 
produced without altering the rotor configuration, for all cases tested. 

' King was further able to use the NPS modeler to simulate the performance of a 
1/6 scale Comanche rotor system. The goal of the simulation was to match the fixed 
system damping values to the test data and to UMARC (computer simulation modeler 
designed at University of Maryland). In this comparison the NPS modeler shows slightly 
better results than UMARC at the lower RPM’s tested. At higher rotor speed, the 


difference between the two simulations is negligible, Figure 6. 
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Figure 6. Final NPS simulation results with increased geometric gain [From Ref. 5] 


From these.results King concluded that “the NPS modeler shows potential in 
modeling rotors with nonlinear physical parameters and further development should be 


considered” [Ref. 5]. 
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Ill. MATLAB® 


A. INTRODUCTION AND HELP 

The NPS modeler has been shown to be effective in modeling coupled rotor- 
fuselage motion. In order to help facilitate the functionality of the simulation program, 
the inclusion of a symbolic processing toolbox in MATLAB provides a means of refining 
the simulation process. The computer programming contained within this thesis is 
accomplished entirely in the computer programs MATLAB® and SIMULINK®. Both 
programs are developed by Mathworks, Inc. and are run interchangeably. MATLAB’ is 
an array and matrix based program allowing programming features similar to other 
computer programming languages. In addition to its matrix orientation, MATLAB® also 
offers Graphical User Interface (GUI) tools and an excellent graphics package. 
MATLAB? stores all data as arrays, which lends itself to the manipulation of sets of data 
in a wide variety of ways. Recently Mathworks, Inc. has incorporated the Symbolic Math 
Toolboxes into MATLAB®’s numeric environment. These toolboxes supplement 


MATLAB®”s numeric and graphical facilities with several other types of mathematical 


computation: 









Facilit 


[| sFacility” 
series 
decomposition, and canonical forms of symbolic matrices 


Simplification Methods of simplifying algebraic expressions 






quations differential equations 
Solution of Symbolic and numerical solutions to algebraic and 
Equations differential equations 


Variable-Precision Numerical evaluation of mathematical expressions to an 
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Arithmetic specified accurac 


transforms 
Functions 
(®) 


Table 1. MATLAB” Symbolic Math Toolboxes [After Ref. 13] 






The computational engine underlying the toolboxes is the kernel of MAPLE®, a 
system developed primarily at the University of Waterloo, Canada, and, more recently, at 
the Eidgendssiche Technische Hochschule, Ziirich, Switzerland. MAPLE is marketed 
and supported by Waterloo Maple, Inc. ) 
Maple V Release 4 was incorporated into MATLAB® 5. There are tow toolboxes. 
The basic Symbolic Math Toolbox is a collection of more than one hundred MATLAB® 
functions that provide access to the Maple kernel using a syntax and style that is a natural 
extension of the MATLAB® lan guage. The basic toolbox also allows you to access 
functions in MAPLE®”’s linear algebra package. The Extended Symbolic Math Toolbox 
augments this functionality to include access to all nongraphics MAPLE® packages, 
MAPLE® programming features, and user defined procedures. With both toolboxes, you 
can write your own M-files to access MAPLE® functions and the MAPLE® workspace 
[Ref. 13]. 
In order to determine it you have the Extended Symbolic Math Toolbox a simple 
check on the MATLAB® command line will work. The syntax is as follows: 
>>maple(with’,’numtheory’) 

This will either: 
1) list a set of libraries that are loaded or 
2) produce the following error: 


222 ’with’(...)’ requires extended symbolic toolbox. 
Error in ==> Paul P.:Applications :MATLAB 4.2a:Toolbox:Symbolic 
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Installer:symbolic:maplemex.mex 

Error in ==> Paul P.:Applications : MATLAB 4.2a:Toolbox:Symbolic 
Installer:symbolic:maple.m 

On line 84 ==> [result,status] = maplemex(statement); 

More information about obtaining the Extended Symbolic Math Toolbox may be 
obtained from the Mathworks website [www.Mathworks.com]. 

This thesis is written assuming the reader possesses a basic understanding of 
MATLAB®. Additional help with the Symbolic Math Toolbox may be obtained in 
several ways. General information about symbolic functions may be obtained by typing, 

>> help sym/function. 

Here, function, is the name of an “overloaded” MATLAB® numeric function. 
This method provides symbolic-specific implementation of the function, using the atte 
function name. This provides some consistency in the logically programming process 
and the functions used in both MATLAB® and MAPLE®. For example diff is a function 
that is used numerically as well as symbolically by MATLAB®. 

>> help diff 
will produce a different result that 

>> help sym/diff 

Help for the numeric version informs you that the function name is “overloaded”. 
Overloaded methods 

>> help char/diff.m 

>> help sym/diff.m 

Help may also be obtained on MAPLE® commands from the MATLAB® work 


space. 
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>> mhelp diff 


This returns the help page for the MAPLE® diff function. 


B. IDENTIFYING SYMBOLIC VARIABLES IN MATLAB® 


Using the Symbolic Math processing capabilities of MATLAB’ is different than 
using MAPLE®. Of particular difference is the requirement that a variable be identified 
as a symbolic variable before it used in MATLAB®. If you have ever received the error 

???Undefined function or variable 7’ 
you have fallen victim to this peculiarity. The exception to this rule is if a variable is the 
result of an equation consisting of previously identified symbolic variable. Then the 
variable is automatically assumed to by symbolic. 

A symbolic variable may be identified in at least two methods. To designate one 
variable as a symbolic expression use 

u=sym(‘u’ 
or 

b = sym(‘beta’). 
Here ‘u’ is established as a symbolic variable u and b as beta. These variables are 
defaulted to be with respect to x. A variable my also be assigned to an expression, 

f = sym(‘a*x42 + b*x +c’) 
but in this form the Symbolic Math Toolbox does not create variables corresponding to 
the terms of the expression, a, b, c, x. To perform symbolic math operations (e.g., 


integration, differentiation, substitution, etc.) on f, each variable must be created 
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explicitly. This may be accomplished in the form above or by typing, 
symabcx. 

Again these variables are defaulted to be with respect to x. The independent variables is 
chosen with the idea that they are typically lower-case letters found at the end of the 
Latin alphabet (e.g., x, y, or z). Therefore the closest letter to ‘x’ alphabetically is chosen 
as the default symbolic variable. If there are two equally close, the letter later in the 
alphabet is chosen. Establishing a variable with respect to another variable may be 
accomplished in the following format, 

zet = sym(‘zet(t)’). 
The variable zet is now defined with respect to ‘t’. This is of particular importance when 
performing integration and differentiation as in this thesis. Further explanations of 
defining symbolic variables in MATLAB® may be found in reference 13 symbolic 


toolbox. 


Cc. SUBSTITUTIONS 

When processing expressions using the Symbolic Math Toolbox, the MAPLE® 
kernel, accessed by MATLAB®, does not “a priori”, treat the expression as a number. 
MAPLE® assumes that the symbolic variables as “a priori” indeterminate. That is, they 
are purely formal variables with no mathematical properties. Consequently, when 
calculating an expression, the variables do not automatically assume a numerical value, 
as assigned. For example, 

>> syms x yZ 


>>l=x+y*z 
>>x=1; 
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>> y = 3; 

>>z=5; 

>>] 
produces the answer, 

l=x*2+y*z 
However, 

1 = subs(l) 
sreaices 

16 
In many cases a specification of the desired form of output is required. If a numerical 
output is desired for the expression and it requires processing beyond basic arithmetic, 
the ‘subs’ function may be incorporated with the ‘double’ command to produce a double 
precision array. Double precision array is the general MATLAB® workspace parameter. 
Inquiry into the class of each variable may be accomplished by typing, 

>> whos. 


A list of the Name, Size, Bytes, and Class of all variable accessible to the workspace are 


listed. 


D. MAPLE® FUNCTION 


The ‘maple function’ allows the user to access MAPLE® functions directly from 
the MATLAB® workspace. This function takes sym objects, strings, and doubles as 
inputs and returns a symbolic object, character string, or double corresponding to the 


class of the input. 
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It is important to know the syntax of the calling sequence used by MAPLE® 
function to be used. This information may be obtained by using ‘mhelp.’ The general 
format for using the MAPLE® function is, 

A(i,j) = maple(‘coeff’,EOM(i),ddDOFq(j)), 
where ‘coeff’ is the maple function to be used. In this case, ‘coeff’ sets A(i,j) equal to the 
ddDOF(j) coefficients of the expression EOM(i). In some instances the format of the 
MAPLE® function is not consistent with form used in MAPLE®. For example the 
MAPLE® function ‘union’ requires the sets to be evaluated to be on either side of the 
function instead of in parentheses after it: 

setl := setA union setB 
instead of, 

setl = union(setA,setB). 

These forms are not conducive to using the maple function however the proper 
form by be accomplished by using strings and converting the syntax. An example of this 


is provided in Appendix I with the MATLAB® function munion. 
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IV. MAPLE® TO MATLAB® CONVERSION 


A brief explanation of the programming tools required to produce the results of 
this thesis is provided above. Unfortunately, excessive trial and error went into 
establishing the differences in the programming format required to produce the desired 
output. The similarities in the two languages proved to be a hindrance rather than 


beneficial. 


A. DEFINING THE VARIABLES 

In order to accomplish the derivation of the equation of motion using Lagrange’s 
equation in MATLAB® the following hurdles need to be crossed. MAPLE® allows 
variables to be sipinden 1.6. 

-W=Qt+¢, 
where @, is subindexed with k, as is commonly used in mathematical text books. In this 


case k represent an individual rotor blade and for the three bladed model k = three. This 
mathematical formulation is represented as 

psi:=Omega*t+Phi[k]; 
in MAPLE®. MATLAB® uses square brackets, [ ], to represent an array or matrix. Asa 
result there is no way to represent k in the MATLAB® workspace as it is used in 
MAPLE®. A multidimensional array was used to resolve this conflict. A third index is 
used to represent the values of each blade. In this form MATLAB® calculates ‘psi’ three 
time for a three bladed model. The syntax in this form is, 


Phi = [Phil Phi2 Phi3]; 
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Psi = Omega*t + Phi; 
But when the variable is called, it is done so with a for loop from | to k (k represents the 
number of blades). 
fori = 1:k 
ml(:,:,1) = [cos(psi(i) sin(psi(i)) 0 
-sin(psi(i)) cos(psi(i)) 0 
00 1); 
end; 
ml is generated 3 times with a different psi each time. This results is psi with Phi, Phi2, 


and Phi3 in each successive generation. This accomplishes the necessary ‘place holding’ 


required for each blade’s energy terms to accomplish the derivation. 


B. SET MANIPULATION 

As a result of MAPLE®”s dependence on taking derivatives with respect to 
independent variables, the time function notation, which represents the degrees of 
freedom and the time rates of change of the degrees of freedom, must be converted to 
independent variable notation as noted in reference 4. To accomplish this a substitution 
of sets was used to switch the necessary terms before and after the differentiation. In 
MATLAB*® this may be accomplished in two ways. The first is similar to the form used 
in the original derivation. 

DOFF gS ees an 

DOFB [zetl, zet2, zet3]; 

DOF = [DOFF, DOFB] ; 


GDOF = diff (DOF,t); 
ddDOF = diff (dDOF,t); 


syms ql q2 q3 q4 q5 dq1 dq2 dq3 dq4 dq5 ddq1 ddq2 ddq3 
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dadgq4 ddq5 


DOFg = [ ql q2 q3 q4 Q5]; 

daDOFg = [ dql dq2 dq3 dg4 dq5]; 
ddDOFg = [ ddql ddq2 ddq3 ddq4 ddq5]; 
setl = []; 

set2 []; 


for i=1:5 % size of DOFg vector 


setl = [setl maple(’‘=‘’,DOF(i),DOFq(i)) J]; 
setl = [setl maple(’‘=‘’,dDOF(i),dDOFq(i)) ]; 
setl = [setl maple(’‘='',ddDOF(i),ddDOFq(i)) ]; 
set2 = [set2 maple(’‘=‘',DOFq(i),DOF(i)) ];. 
set2 = [set2 maple(’‘=*‘',d@DOFq(i),dDOF(i)) ]; 
set2 = [set2 maple(’*=‘’,ddDOFq(i),ddDOF(i)) ]; 

end 

setl = maple(’convert’,setl,’set’); 


set2 = maple(’convert’,set2,’set’); 

This method defines the terms required as symbolic variables. It defines the sets 
and then uses a maple function to set the terms equal to each other for later substitution. 
It then uses another maple function to convert the sets into a form that is recognizable by 
MATLAB®. 

The second method uses the ‘subs’ command local to the MATLAB® workspace. 
By establishing the necessary variables into an array of ‘old’ and ‘new’ terms they can be 
switched term for term in a given variable (array). 

DOFF = [ul,u2]; 

DOFB = [zetl, zet2, zet3]; 

DOF = [DOFF, DOFB]; 

GDOF = diff(DOF,t); 

daDOF = diff (dDOF,t); 

tempDOF = [DOF dDOF ddDOF] ; 


syms ql q2 q3 q4 q5 dql dq2 dq3 dq4 dq5 ddgql ddq2 ddq3 
ddq4 ddq5 


DOFg = [ ql q2 q3 q4 q5]; 
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aDOFg = [ dql dq2 dq3 dq4 dq5]; 

ddDOFg = [ ddql ddq2 ddq3 ddq4 ddq5]; 

tempDOFq = [DOFq dDOFq ddDOFq] ; 

Temp = subs(T,dDOF,dDOFq) ; 
T, ul, u2, zetl, zet2, zet3 are predefined arrays and tempDOFq is used in a later portion 
of the derivation. The speed required to run each method varies depending on the 
operation to be accomplished. For a simple substitution, the MATLAB® ‘subs’ function 
proved to be the quicker method but when larger symbolic terms were involved the first 
method proved to be the quicker solution. Therefore an optimization was done to 


determine the optimum combination of the two methods in order to minimize the 


computational time required to complete the simulation. 


C. STATE DERIVATIVE CALCULATION 


Using the Symbolic Math Toolbox in MATLAB allows calculation of the state 
derivatives without first creating a ‘B’ array, as is necessary with the MAPLE® ersien: 
The B array is used to convert the output to C or Fortran code, before it is placed in a 
MATLAB? function. By avoiding this step, both the MATLAB?® function and the 
SIMULINK® faite can cal] the state derivative directly from the workspace. This is 
extremely beneficial in that it eliminates any errors associated with converting generated 
code to a format that usable in SIMULINK®.. It eliminates extra steps in the setup 
process, helps facilitate simulations, and reduces the time taken to convert the generated 


code. 


34 








D. VARIABLE ROTOR BLADE MODELING 


After the code was found to run successfully in SIMULINK®, it was further 
adjusted to allow for a variation in the number of rotor blades. In the MAPLE® version 
this was self-contained by using the subindex ‘k’. In the MATLAB® version, without the 
capability of subindexing, each array was built up as a temporary array, capable of 
deriving a seven bladed model, then only the necessary terms were used to complete the 


derivation. 


E. COMPLEX MODEL 

The complex model as derived by Robinson, [Ref. 4], was converted into a 
MATLAB® function, Appendix E. This model is very applicable to future helicopter 
design due . the added degrees of freedom. It is the complex model and variations of it 
that will allow the designer to examine the stability characteristics of rotor systems 
without expensive and time consuming tests in a wind tunnel. The flexibility in alteing 
the degrees of freedom and accessing the equations of motion with each time step are 
what make the simulation useful beyond the calculations performed using the simple 
model. Additionally the inclusion of nonlinear terms opens a whole new field of 
dynamic modeling to be considered. Some of the possible areas of exploration are 
damper springs, duffing springs, and the stability characteristics of varying the length 


rotor blades. 


F. MAPLE® PROCEDURE 


Reference 13 describes a method to input a MAPLE® symbolic derivation into 


the MATLAB® workspace. Creating a source file using a ‘maple procedure’ the 
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Extended Symbolic Math Toolbox is capable of reading the source file. Unfortunately 
the results produced little to no gain due to the process of changing the source code 
identifier. This process was further complicated due to windows automatically 
associating the file as a text file or attaching its own identifier to file. Accessing the 
sources code was eventually possible by changing the file in DOS mode. This process 
could possibly be simplified using a UNIX based system however it does not produce the 


ultimately desired results of this thesis. 


36 








V. SIMULATION MODEL 


A. ORIGINAL SIMULINK® S-FUNCTION 

Introduction of the rotor-fuselage dynamics into the Simulink® simulation is 
accomplished through the use of a Simulink® S-function. The S-function is a user 
defined Simulink® block that outputs the rotor time histories as functions of the user 
defined rotor and fuselage parameters. SIMULINK® is a convenient Graphical User 
Interface (GUD) that performs simulation time integration, allows user input changes 
between simulations, and provides the user with graphical output. 

SIMULINK® accesses an S-function through is numerical integration routines. 
The routines make calls to the S-function for specific information, the type of information 
returned is dependent on the value of a flag variable sent by the cuepmianl routine. For 
example, 

flag =O S-function returns sizes of parameters and initial conditions, 

flag = 1 S-function returns state derivatives dx/dt, 

flag =3 S-function returns outputs. 
The section of the S-function which computes the derivatives at each time step is the 
section which calls the equations of motion. 

The Lagrangian derivation of the equations of motion generated the equations in 


the form, 
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where X is a vector of displacement degrees of freedom of the system. The output of the 


derivation was further manipulated into the equivalent form, 


l4a@z.oF = 7G.z.0 
where A is an NxN matrix and f is an Nx1 vector, with N = number of degrees of 


freedom of the system. This is possible since the equations are quasi-linear in the second 


derivative terms, i.e., no terms of types such as X, or sin( x ), etc. This form can then be 
transformed from N second order equations to 2N first order equations as follows, 
x=w 

w=[4P’f 

These equations are evaluated at each time step in a numerical simulation to give 
the state derivatives. This study uses a Runge-Kutta algorithm in SIMULINK® to solve 
the numerical ordinary differential equations. The Runge-Kutta algorithm generally 
outperforms other schemes for systems of nonlinear ordinary differential equations which 
are not too stiff and they handle discontinuities well [Ref. 14]. 

A Coleman-like 3-bladed hub-rotor model was used for the initial simulation 
model. It provided for linear stiffness and damping, quadratic damping, and cubic 
stiffness in the blade degrees of freedom. Other, simulations for damperless rotor 
analysis have been created for four-blade and five-blade hub-rotor simulations. Table 2 
provides a break down of the variables used in the S-function and Figure 7 shows the 


SIMULINK® model utilized for the simple mode] simulation. 
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Table 2. Parameter inputs for simple model 





Fuselage Displacement 


> 3-blatie rotor-pyton (NL) Demux 


Helo 


Mux — | 


Mux Lead-lag Displacement 


Demux 


Figure 7. SIMULINK® diagram of NPS simple model 
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B. MODIFICATIONS TO SIMULINK® S-FUNCTION 


Of significant importance in this and future versions of the simulation is 
MATLAB*”s definition of variables. From the original MATLAB® function where the B 
matrix was input as optimized C or Fortran code, the variables in the equations of motion 
were in array format. Using the Symbolic Math Toolbox requires the variables to be 
defined as symbolic. In doing so, an additional substitution step is required to solve the 
state derivatives as real variables. This substitution may be accomplished in either the 
derivation function or the MATLAB® S-function. A coordination between the two 
functions is required in order to properly derive the equations of motion and still run the 
simulation. For example time (t) is needed to derive the equations of motion, however it 
is not recognized in MATLAB® unless it is defined as a symbolic variable prior to the 
derivation, but once it is defined as a symbolic variable an actual value for ‘t’ (each time 
step) may not be substituted in to the simulation model until it is redefined as a ‘double 
array’ variable. 

Using a MATLAB® function to derive the equations of motion allows the model 
to call the derivation function with each time step. Unfortunately this process is very 
time consuming and severally slows the simulation process. The potential exists to cut 
and past the equations of motion into an S-function, in order to reduce the simulation run 
time. This method was not used in an attempt to maintain the integrity of the function 
calling process. Future thesis may have the opportunity to optimize this process. It 
should be noted that the cut and paste method would still significantly reduce the time, 


labor, and possible errors associated with converting the MAPLE® code. 
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As an alternative mere deriving the equation of motion with each time step, a 
data file was created with the state derivatives input from the derivation function. The 
goal of this setup was to reduce the computational time. This method was not found to be 
as effective as anticipated because the process still required the derived variables to be 
converted to “double array’ variables before value substitution in the S-function 


workspace. 


C. USING THE SIMULINK MODEL 

Six files are required to run the simulation. For the three bladed simple model the 
files are: simple3.m, helo3b5.m, simple.m, simple.mdl, signum.m, signum1.m, abs.m. 
From the MATLAB® command window type in 

>> simple 
a figure will appear identical to Figure 7. The operator may begin the simulation by 
“clicking” on the arrow @ ) on the menu bar. Two figures will be displayed, Fuselage 
Bieieenent and Lead-lag Displacement. The Fuselage Displacement figure represents 
the position of the fuselage in the x/y plane with respect time. The Lead-lag 
Displacement figure represents each blades position with respect to time. From here the 
operator can deduce whether the model is stable or unstable. 

The rotor parameters may be changed by “double-clicking” on the box titled “3- 
blade rotor-pylon {NL}. A new window will appear containing the “Block Parameters.” 
Each of these terms may be changed to effect the desired response from the model. An 


explanation of these terms is included in Table 2. Further explanation of using the NPS 


simulation model may be obtained in reference 17. 
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VI. RESULTS 


A. EQUATIONS OF MOTION 

The equations of motion calculated for the simple model by MATLAB®, 
Appendix A, are equivalent to the equations calculated by Robinson [Ref. 4]. As a result | 
the plots produced by SIMULINK® display the same time history plots, Figures 8, 9. 
The derivation of the equations of motion using MATLAB® Symbolic Math Toolbox are 


complete. 


Hub Motion in Horizontal Plane 





Hub y-axis displacement (ft) 
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Figure 8. Fuselage displacement for unstable region 
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Lead-Lag Displacement Time Histories 
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Figure 9. Lead-Lag displacement for unstable region 
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VII. CONCLUSIONS AND RECOMMENDATIONS 


Previous work by Robinson and King was based on a derivation in which the 
equations of motion were derived in MAPLE® then converted into Fortran or C where the 
equations of iiss were further converted in a MATLAB® usable format (Method A). 
The focus of this thesis was to derive the equations of motion using the Symbolic Math 
Toolbox in MATLAB®(Method B). This derivation allows direct communication 
between the derivation function and the simulation thus eliminating user interface in the 
simulation process. 

Based on the simulation process utilized in this thesis, Method B, deriving the 
equations of motion using MATLAB® was found to possess both advantages and 


disadvantages. 


ADVANTAGES 
- No required conversion of equations of motion to MATLAB® readable code 
- Eliminates “user” interface 


- Single software package 


DISADVANTAGES 

- Time consuming 

- Symbolic processing format is not as usable as MAPLE® 
- Does not provide instant feedback to equation processing 


- Requires additional substitution step 
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TRADEOFF 

- It was found, using Method B, that eliminating user interface from the simulation 
process incurs a substantial run time penalty 

- With minor user interface, modification to Method B eliminates the time penalty 
incurred 

- The ino user interface still reduces the potential for error associated with code 
conversion 


- Modification to Method B maintains a comparable run to Method A 


Expanding on the “Tradeoff Option”, if the user intends to utilize the NPS modeler in the 
form in which the simulation calls the derivation function internally, Method B, then a 
run time penalty is incurred. However if modifications to this process are acceptable, the 
result is a reduction in potential errors associated with code conversion and a savings in 
the time required to perform this process. 

For example, potential modification to the simulation process used in Method B is 
the option of “cutting” the resultant matrices from the MATLAB® workspace and 
“pasting” them into the simulation workspace. This “cut and paste” method is similar to 
deriving the equations of motion in Method A. However, it eliminates the process of 
converting the code from Fortran or C to MATLAB?® executable code thus reducing 


“user” interaction and eliminating the potential for user contamination of data. 
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Future emphasis should be placed on optimizing the passing of variables between 
the MATLAB® derivation function and the simulation process. Continued analysis of 
this process may determine a more effective means of substituting the required variables 
and reduce the time to run the simulation. With further optimization of the NPS modeler 
in the MATLAB® format and the continued increase in computational processing speeds, 


eliminating user interaction between the derivation and simulation process is a realistic 


goal. 
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APPENDIX A. MATLAB® EQUATIONS OF MOTION FOR SIMPLE INPLANE 
(COLEMAN) MODEL WITH THREE ROTOR BLADES 


pretty(EOM) 


[ ps 
[1/2 v1 dql abs(1, dql) + vl dq1 | dql |+M1 ddqi 


2 
- 2mb2 R %4 Omega cos(q4) dq4 - mb2 R %4 Omega cos(q4) 


2 2 
-mb1 R %6 Omega cos(q3) - mb1 R %6 cos(q3) dq3 


2 
- mb1 R %6 sin(q3) ddq3 + mb2 R %3 Omega sin(q4) 


2 
+ mb2 R %3 sin(q4) dq4 - mb2 R %3 cos(q4) ddq4 


2 2 
- mb3 R %2 Omega cos(qg5) - mb3 R %2 cos(q5) dq5 


2 
- mb3 R %2 sin(q5) ddg5 + mb1 R %5 Omega sin(q3) 


2 
+mbl1 R %5 sin(q3) dq3 - mb3 R %1 cos(q5) ddq5 


2 
- mb1 R %5 cos(q3) ddq3 + mb3 R %1 Omega sin(q5) 


Z 2 
+mb3 R %1 sin(q5) dq5 - mb2 R %4 cos(q4) dq4 


- mb2 R %4 sin(q4) ddq4 + K1 ql + mb1 ddql + mb2 ddq1 +1 dql 
+mb3 ddq1l + 2 mb1 R %5 Omega sin(q3) dq3 
- 2 mb1 R %6 Omega cos(q3) dq3 + 2 mb3 R %1 Omega sin(q5) dq5 


~ 2 mb3 R %2 Omega cos(q5) dq5 + 2 mb2 R %3 Omega sin(q4) dq4 
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2 2 Z 
- mb3 el %2 Omega - mb2 el %4 Omega - mb1 el %6 Omega , 


2 2 
1/2 v2 dq2 abs(1, dq2) + v2 dq2 | dq2 | - mb3 R %1 Omega cos(q5) 


2 
-mb1 R %6 sin(q3) dq3 + mb1 R %6 cos(q3) ddq3 


2 
-mb1 R %5 Omega cos(q3) - mb1 R %5 sin(q3) ddq3 


2 2 
-mb2 R %3 Omega cos(q4) - mb1 R %5 cos(q3) dq3 


+ mb2 R %4 cos(q4) ddq4 + mb3 R %2 cos(q5) ddq5 


2 2 
- mb3 R %2 sin(q5) dq5 - mb3 R %2 Omega sin(q5) 


2 
-mb3 R %1 cos(q5) dg5 - mb2 R %3 sin(q4) ddq4 


2 2 
-mb1i R %6 Omega sin(q3) + mb1 ddq2 - mb2 R %3 cos(q4) dq4 


2 
- mb3 R %1 sin(q5) ddg5 - mb2 R %4 sin(q4) dq4 


2 2 Z 
-mb2 R %4 Omega sin(q4) - mbl el %5 Omega - mb3 el %1 Omega 


+ mb3 ddq2 + mb2 ddq2 - 2 mb3 R %1 Omega cos(q5) dq5 

- 2 mb1 R %5 Omega cos(q3) dq3 - 2 mb1 R %6 Omega sin(q3) dq3 
- 2 mb3 R %2 Omega sin(q5) dq5 - 2 mb2 R %4 Omega sin(q4) dq4 - 
- 2 mb2 R %3 Omega cos(q4) dq4 + M2 ddq2 + K2 q2 + c2 dq2 


2 2 
-mb2 el %3 Omega , mbl el Omega R sin(q3) + Czetal dq3 + Kel q3 


54 


% 


— 


%2 


%3 


GA 


%5 


%o6 





3 2 2 
+Kd1 q3 + Vzetal dq3 abs(1, dq3) - ul(t)+mb1R ddq3 


+ mb1 ddq2 R %6 cos(q3) - mb1 ddq2 R %5 sin(q3) 
- mb1 ddq1 R %6 sin(q3) - mb1 ddq1 R %5 cos(q3) 


2 
+2 Vzetal dq3 | dq3 |, Vzeta2 dq4 abs(1, dq4) 


3 
+2 Vzeta2 dq4 | dq4|+Kd2q4 +Ke2 g4 - u2(t) 


2 
+mb2 el Omega R sin(q4) - mb2 ddq1 R %4 sin(q4) 


2 
- mb2 ddq1 R %3 cos(q4) + mb2 R ddq4 + mb2 ddq2 R %4 cos(q4) 


| 
- mb2 ddq2 R %3 sin(q4) + Czeta2 dq4 , Czeta3 dq5 + Ke3 q5 + Kd3 q5 


” 
- u3 + Vzeta3 dq5 abs(1, dq5) + 2 Vzeta3 dq5 | dq5 | 


2 
- mb3 ddq1 R %2 sin(q5) - mb3 ddqi R %1 cos(q5) + mb3 R ddq5 


+ mb3 ddq2 R %2 cos(q5) - mb3 ddq2 R %1 sin(q5) 


5 
+mb3 e1 Omega R sin(q5)] 


== sin(Omega t + Phi3) 
= cos(Omega t + Phi3) 
= sin(Omega t + Phi2) 
== cos(Omega t + Phi2) 
:= sin(Omega t + Phil) 


:= cos(Omega t + Phil) 


a5 
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APPENDIX B. MATLAB® WORKSHEET FOR SIMPLE INPLANE (COLEMAN) 
MODEL WITH THREE ROTOR BLADES 


= simple3 () 


£1) 


[Al, 


function 


EQUATIONS OF MOTION FOR A HELICOPTER IN GROUND RESONANCE 


CONSIDERING ONLY INPLANE DEGREES OF FREEDOM 


Xe) 


re) 


° 
6 


clear 
cle 


bles for transformations 


varia 


Definition of 


syms Phil Phi2 Phi3 Omega t 


[Phil Phi2 Phi3]; 


Phi 
psi 


. 
, 


a: 


Omega*t + Ph 


sym('zetl(t)'); 


zetl 


zet2 = sym(’zet2(t)’); 


zet3 = sym(’zet3(t)’); 


. 
‘ 


[zetl zet2 zet3] 


zet 


% COORDINATE TRANSFORMATIONS 


0 


[cos(psi(i)) sin(psi(i)) 


(psi (i)) 


0 


cos(psi(i)) 


-sin 


00 1]; 


[cos(zet(i)) sin(zet(i)) 0 


:,2,1) 
~sin 


m2 ( 


cos(zet(i)) 0 


(zet (i)) 


~ 


= mi(:).’ 


=m2(:).'; 


m2 (:) 


oe 


ENERGY OF ROTOR BLADES 


oe 


de 
oe 


Kinetic energy of rotor blade (TBk) 


ce) 


Xx 
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syms el R mb1l mb2 mb3 


’ 


[mb1 mb2 mb3] 


ul = sym(’ul(t)’); 


‘ 


[ul u2 0] 
fel 0 0] 
[R 0 0] 


rhoHI_I 
rhoBuH 
rhoPBd = 


, 


. 
a 


rare 


7 


rhoBuH*ml ( 


rhoBuH_I 
rhoPBd_I 


vg Ls 


ia 


Sok) m2 ¢ 


. 
Cae i 


rhoPBd*mi ( 


a 
tH 
Ha 
O 
a 
re 
+ 
a 
ye) 
rea) 
AY 
fe) 
‘G 
re 
+o 
z . 
mp 
3 s 
fo oO 
OG 
GH 
Sd 
~ WY 
Su 
ed 
do OU 
vy il 
& 
od An 
Y -d 
How 
pe 
“> 


re On eee 


= 1/2*mb(i)*(V(1,:,1)%2 + V(2, 


2,1) 


esr 


TBk ( 


end 


(UBk) 


% Potential energy or rotor blade 


syms Kel Ke2 Ke3 Kdl Kd2 Kd3 Ksl Ks2 Ks3 z 


[Kel Ke2 Ke3]; 


[Kd1l Kd2 Kd3]; 


[Ks1 Ks2 Ks3]; 


Ks 


for i 


1/2*Ke (2) *zet (1) 2 


a 


UBk1 


4 


1/4*Kd(i)*zet(i)%*4 
UBk3 = 1/4*Ks(i)*signum(zet(i)-z)*(zet(i)*%2+z%*2- 


2*zet (i) *z)+1/4*Ks (i) *signum(zet (i)+z) * (-zet(i)*2-z%*2- 


2*7eb (i) *2)41/2*Ks (1) *zet (2) °2. + 1/2*Ks(i) *272 


UBk2 


, 


UBk1 + UBk2; 


UBk (i) 


end 


(DBk) 


Dissapative function of rotor blades 


° 
© 


syms DB1 DB2 DB3 Czetal Czeta2 Czeta3 Vzetal Vzeta2 Vzeta3 


. 
a 


[Czetal Czeta2 Czeta3] 
[Vzetal Vzeta2 Vzeta3] 


Czeta 
Vzeta 


cA 
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Lia Czeta (2) (dal fh(2eu(1) 7) 72:4 


) 


i 


DBk ( 


)*(difi(zet(i),t))*2*abs (aiff (zet(i),t)); 


i 


Vzeta ( 
end 





oe 
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oe 


HUB 


ENERGY OF 


oe 


ae 
oe 


/ Dissapative 


Potential energy (UH) 


iy 


Kinetic energy TH) 


function 


° 
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(DH) 


syms Kl K2 cl c2 vl v2 M1 M2 
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, 


[cl c2] 


aang 
Pam 
ddd 
aay 
iow 

Pam 


. 
, 


EECUZ (6772 


an 


L/2°M (1) dt FEU, b)o2 + b/2*M (2) 4d 
Loko) Paes. OL se 2 Ro) Pe 2 


DE ‘= 1/2 etl) Carer Gil, by yes 


TF 


. 
f 


UF 


1/2*e(2)* (diff (u2,t))*24+1/2*v(1) * (diff(ul,t))*2*abs (diff (ul 


pel EL 2* (2) * (ALPE (a2 2b) ) 2" abs (dati tu2,c))> 


d displacements 


LZe 


’ 


Generalized forces on general 


% 


syms Fl F2 F3 F4 F5 u3 


{ful u2 u3]; 


ust 


Fl 
F 


[PL PZ. ES 2a Pods 


and zet(2) 


; 6 or zet(l1) 


zet (3) ] 


zet(2), 


= [zet(1), 
[DOFF, 
ADOF = diff (DOF,t) 


DOFB 
DOF 


° 
cf 


DOFB] 
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cA 


diff (dDOF, t) 


ddDOF = 
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, 


te 
] 


] 
] 


1? 
if 


a 





. 
, 


ta 


‘* ,AadDOF (i) , dADOFgq(i) ) 


‘* , ,DOFq(i) ,DOF (i) ) 


a 


ow 


) 


ies 


size of DOFq vector 
{[setl maple(’‘=*‘’,DOF(i),DOFq(i)) 


[setl maple(’*=*‘'’,dDOF(i),dDOFq(i) ) 
[set2 maple(’ ‘=*‘’,dDOFq(i),dDOF(i) ) 
[set2 maple(’ ‘=*’,dd@DOFq(i),ddDOF (i) ) 


[set2 maple(’’ 


[ ddql ddg2 ddq3 ddq4 ddq5] 
[setl maple(’* 


{ A@gql dq2 dq3 dq4 dq5] 


maple(’convert’,setl1,’set’) 


[ ql q2 q3 q4 q5] 
set2 = maple(’convert’,set2,’set’); 
= size(DOFB) ; 


setl 
setl 
setl 
set2 
set2 
TF; 
UF; 
= DF; 


set2 


syms ql q2 q3 q4 q5 dql dq2 dq3 dq4 dq5 ddql ddq2 ddq3 ddq4 
Substition set construction 


ddq5 
DOFq = 
aDOFg = 
ddDOFq = 
end 

setl 

D1 

[ln] 


U = U+UBk ( 
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a 


f 


D1+DBk (i) 
maple(‘subs’,setl1,T) 


= size(DOF); 


D1 
Switch of time dependent terms 


Complete EOM 


end 
Temp 
[1,n] 





‘diff (Temp, ADOFq(i)) 





. 
i 


templ 
temp2 
temp3 


f 


maple(’subs’,set2,temp1) 


diff (temp2,t) 
maple(’subs’,set1,temp3) 


diff (Temp, DOFq(i)); 


tempL3 


‘ 


° 
f 


L2 


maple(’subs’,seti,U); 


check dimensions of DOFq 


maple(’subs’,set1,D1) 


ae 
Keo} 


diff (tempL3,DOFq(i)); 


tempL4 


L3 


. 
f 


check dimension of DOFgd 


% 


. 
oa 


simplify (L1-L2+L3+L4-F(i)) 


diff (tempL4, dDOFq(i) ) 


EOM (i) 
of EOM and F 


end 


L4 


check dimension 


% 


. 
i 


Elimination of excess terms 


% 


[ signum1(1,q3-z) signuml(1,q4-z) signuml1(1,q5-z) 


signuml1 (1,q3+z) signum1(1,q4+z) signuml(1,q5+z) 


temps 


maple(’abs’,1,dq3) maple(’abs’,1,dq4) maple(’abs’,1,dq5) 


i 


maple(’abs’,1,dq2) 


L- 0-0-0 0. 0.0: 0° OO 0-0] 


maple(‘abs’,1,dql1) 


temp0 = 


. 
, 


‘’ ,tempS (i) ,temp0(i)) 


{setS maple(’‘= 


sets 


end 


. 
a 


maple(‘convert’,setS,‘set’) 


sets 


RATE Al AND £1 FROM EOM 


3% 
GENE 


do odo 


de 


2. 
ie) 


Determine second derivative terms 


maple(’coeff’ ,EOM({(i),ddDOFq(j)); 


maple(’subs’,setS,A(i,j)) 


f 


= ddDOFg*A; 


Ax2dot 


Eliminate second derivative terms from EOM 


fo) 
ie) 
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. 
a 


(EOM (i) -Ax2dot (i) ) 
maple(‘subs’,setsS, f ( 


bd ee 
f£ 


)); 


a. 


i) 


( 


dd); 


lify(f£(i 


imp 


-(s 


end 


Generate state vector 


% 


syms xl x2 x3 x4 x5 x6 x7 x8 x9 x10 


[xl x2 x3 x4 x5]; 


X1ldot 


[x6 x7 x8 x9 x10]; 
[Xltemp Xldot] 


Xltemp = 


X1 


a 


yd 


[DOFq G@DOFq] 


DOFgtemp = 


% Convert terms to state vector form 


. 
f 


] 


‘* ,DOFqtemp (i) ,X1(i)) 


[setX maple(’ ‘= 


setx 


end 


. 
, 


maple(’convert’,setX,’set’) 


setx 


maple(’subs’,setX,A); 


Al 


maple(‘subs’,setX,f£); 


>a 


from input ‘u’ 


fe ta 


Eliminate variable 


. 
‘ 


[ul u2 u3] 
subs (A1,u,ut) 


syms ul u2 u3 


ut 
Al 


f 


. 
i 


subs (f£1,u,ut) 


pe 
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APPENDIX C. SIMULINK® S-FUNCTION FOR SIMPLE INPLANE 
(COLEMAN) MODEL WITH THREE ROTOR BLADES 





function [sys, x0] = helo3b5(t,x,u,flag,I1,12,1I3,14,15,16); 


function [sys, x0] = 
elo3bA(t,x,u, flag,11,12,13,14,15,16) 


S-function arguments: 


t = time 

x = state vector 

u = input vector 

flag = Switch used by numerical integration 


simulation) 
routine to access certain parts of the s- 
unction 


S-function input parameters: 


3 = [mb(1),mbo(2),mb(3),M(1),M(2)] 
I2 =  [R,Omega,el,z] 

I3 = [Phi(1),Phi(2),Phi(3)] 

r4 = [e(1),¢e(2),v{1),v(2), 


Czeta(1),Czeta(2),Czeta(3), 
Vzeta(1),Vzeta(2),Vzeta(3)] 


I5 = [Ketl) Ke 2) Kets); 
Kd(1),Kd(2),Kd(3), 
Ks(1),Ks(2),Ks(3).,. 
K(1),K(2)] 

T6 = (xrXi,xrYi,xrli,xr2i,xr3i, 


xX1i,xXYi,x1li,x2i,x3i] 


dP dP dP DP GP GP GP oP oP oP GP dP dP GP Oo oO dO cP dP dP GO DO DO Fh ach —~ oO oD oO dP dO oP oP DY oo 


S-function to represent dynamics of 3 bladed coupled 
otor- 
fuselage model which considers only inplane degrees of 


oh 
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% freedom, i.e., x and y translational fuselage degrees of 
freedom 
% and lead-lag rotor blade degrees of freedom. 


% Explaination of variables: 





Rn aaa i a eS, eee ORES ay Po 
% 

% mb -> mass of blade 

3M -> effective mass of fuselage 

% R -> distance from lead-lag hinge to blade center 
of mass 

% el -> blade hinge offset 

% Omega -> rotor speed 

% Z -> angle at which blade hits stops 

% Phi -> blade phase angle w.r.t. azimuth postion 

Sc -> fuselage linear damping 

SV -> fuselage hydraulic damping 

%® Czeta -> blade linear damping 

%® Vzeta -> blade hydraulic damping 

% K -> effective stiffness of fuselage (landing 
gear stiffness) 

% Ke -> blade elastic spring constant 

% Kd -> blade duffing spring constant 

% Ks -> blade stop effective spring constant 

S xr_ _i -> initial rate 

6 xX __i -> initial displacement 

% 
LEESESEEEEEESESEESEEEESEEEEEESSEEEEEESESESESEESEEEEESSESEESESSS 
% Define input parameters 
EESEEEESESEESESEESESEETESEESEEEEESEEESEEESSESESESSEEECESEEESESESESS 


S$mbo1=1I1(1) ;mb2=11(2) ;mb3= 
$R=12 (1) ;Omega=12(2);el1=I 
S$Phil=13 (1) ;Phi2=13(2);Ph 
$c1=14 (1); c2=14(2) ;v1l=14(3) ;v2=14(4); 
SCzetal=14 (5) ;Czeta2=14(6) ;Czeta3=14 (7); 
SVzetal=14(8) ;Vzeta2=14 (3) ;Vzeta3=14(10); 
$Kel=15 (1) ;Ke2=15 (2) ;Ke3=15(3); 

SKG1=15 (4) ;KA2=15 (5) ;Kd3=15 (6); 

SKs1=15 (7) ;Ks2=15(8) ;Ks3=15 (9); 
$K1=15 (10) ;K2=15(11) ; 


xrXi=16 (1) ;xrYi=1I6 (2) ;xr1i=16 (3) ;xr2i=16 (4) ;xr3i=16(5); 
xXi=16 (6) ;xYi=1I6(7) ;x1i=1I6 (8) ;x2i=16 (9) ;x3i=1I6(10); 
ESEEECEEEEEEEEEEEEEESEEEEEESEESEEESEEEESSSEEEESEEEEESESEEEES 
% S-function flag conditionals 








1f£ flag == 


sys=[10,0,10,3,0,0]; 
x0=[xrXi,xrYi,xrli,xr2i,xr3i,xXi,xYi,xli,x2i,x3i]; 


elseif abs(flag) == 1 
% Calculate derivatives 


[Al, £1] = simple3; 


xl = x(1); x2 = x(2); x3 = x(3); x4 = x(4); x5 = x(5); x6 = 
x(6); x7 = x(7); x8 = x(8); x9 = x(9); x10 = x(10); 
ul = u(1); u2 = u(2); u3 = u(3); 


x = subs(x); 
u = subs(u); 
fl = subs(f1); 
Al = subs(Al1); 


sys = zeros(1,2*5); 
sys(1:5) = Al\f1.’ 
sys(6:10) = x(1:5); 
% Output states 


elseif abs(flag) == 3 


oe 


psi=Omega*t; 

psil=psi+Phi(1); 
psi2=psi+Phi(2); 
psi3=psi+Phi(3); 
xcl=el*cos(psil)+R*cos(psil+x(8)); 
xe2=el*cos(psi2)+R*cos (psi2+x(9)); 
xc3=el*cos(psi3)+R*cos (psi3+x(10)); 
yel=el*sin(psil)+R*sin(psil+x(8)); 
MS 
) 


oe oP 


oe 


, 


of oO 


oP oP 


oe 


yo2=el*sin(psi2)+R*sin(psi2+x (9) 


de 


yco3=el*sin(psi3)+R*sin(psi3+x(10)); 
xce=(mb(1)*xcl+mb (2) *xc2+mb(3)*xc3)/(s ee 
ye=(mb(1)*ycl+mb (2) *yc2+mb(3)*yc3)/(s ; 
Mag=sqrt (xc*2+yc%2); 
Lf (xe-<.) & yo < 0). | tye = 0) 

theta=pi+atan (yc/xc); 


AP dP of oP 


de 
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% else 


% theta=atan(yc/xc); 
% end 
% gamma=theta+pi; 
% alphal=(psil/ (2*pi)-floor(psil/(2*pi)))*2*pi; 
z alpha2=(psi2/(2*pi)-floor(psi2/(2*pi)))*2*pi; 
% alpha3=(psi3/(2*pi)-floor(psi3/(2*pi)))*2*pi; 
sys (1:10) =x; 
% sys (11)=R*Mag*sin(gamma-alphal) ; 
% sys (12) =R*Mag*sin(gamma-alpha2) ; 
% sys (13)=R*Mag*sin(gamma-alpha3) ; 
% sys (14) =theta; 
else 
sys = []; 
end 
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le() 


simp 


MOTION FOR A HELICOPTER IN GROUND RESONANCE 


WITH VARIABLE NUMBERS OF BLADES 


£1) 


fa 


APPENDIX D. MATLAB® WORKSHEET FOR SIMPLE INPLANE MOTION 


CONSIDERING ONLY INPLANE DEGREES OF FREEDOM 


function [Al, 
EQUATIONS O 


fe) 
6 
% 





4 
10) 
4) 
; ; 9 
vo e) q 
“A “ 
0) oS) O 
O n q vp 
© oO © S 
cd Q, i © 
Q = p 
G n 
fe) n wv 0) qd 
41 Q 5 oO @) 
og © yw Uv 
0) QW -d r ag 
@)) WM oN 1c) COD 

OG © oO 9) eae ee @ 

O-d n G 3 WM WM - 

© G bo: “cl WH g@aqay 

ri HAY Od, ie)) 00 Q 

0 Oo aG4:+Gs aw OU YU 

34 Oo BOA? ono 

rd . 
Wood OBE a & W aae 
O ae) rt} G UH O N tcl ered ert 
wi oO Vv HoT H AT YO “up 
Oo 0 7) Qo od § q aA &% VU 
a u4 OH 3 @ OH nun oO 

n 44 ae 8 CCH YH WH 
On & fe) 0 eo OY ed U0 OW 
Oo g O “of Godoy sv “di G Oo 
© & UY 00 GOH NW GN P-d 
qo yy Od BsNAGwWH M U4 OQ, 
Q 0 gG oO © gO Oo Ow O 

> 0 AY GO OH DP How 
Wea Oo GnooeAaAOOd Gea Vd 
owgd S aD) 

LW ON WDOnrHA WWD oo OW 
ano wv CO OnHnT HV YUOVDTV YW OUT 
QW Nn GH AGNN GC GH Cw W 
CO We HAOfA BSB BAHAY dad 
=O MO) AYoeQHyy QQ Ww QQ Q 

n 

on 
A A A NAAAANAAAAGAAA 
1 ott PopobobPoteotoveoto ot wore ot 

Ww 

“d 

n © GO O Aa) 

a) o pp 10) 

© 0) -d Vv ov 

2 Rad & G NON 4uod n 

aA Lae VOONM VPUP™M i MMM 

dP dP do dP O dP GP GP dO dO de cP DP GO Md GO BD CP 





f 


number of blades 
% blade position angle 


% 





sym(’zet2(t)’) 
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zet2 


‘ 


f 


? 


:k) 


[Phil Phi2 Phi3 Phi4 Phi5 Phi6é Phi7] 


Phit (1 
sym(*zetl(t)*); 


Omega*t + Phi 
sym(’zet3(t)’) 


clear 

cle 

% Define transformation variables 
Phit 

Phi 

psi 

zetl 





zet6 = 


. 
f 


zet5 = sym(’zet5(t)’); 


, 


. 
, 


sym(’zet4(t)’); 

sym(‘zet7(t)’) 

[zetl zet2 zet3 zet4 zet5 zet6 zet7] 
:k) 


zett(l1 


sym(’zet6(t)'); 
zet7 
zett 


zeta 
zet 





de oo de 

oe ct G oe oe 

oe re) (@) oe de 

oe “d -d oe de 

oe 7S) pV 3° de 

ae © © oe de 

de & & oe de 

oe uM uy 3° de 

de e) (@) oe oe 

foe 44 YW oe de 

oe Y 0) oe roe) 

oe G G oe de 

de ~O © de oe 

oe M 4 ae de 

oe zs) :e) oe de 

fond oe oe 

oe ae) ue) oe WY de 

oe (o) YO G oe fx} do 

oe M e) oe A oe 

oe ~~ “do-s 6) fo\e) fd 

oe ~~ Ws oO oe -J dO 

oe “cd ‘d N oe M oo ~ 

de = co oe de 

oe “cd v oe of % oe OM 

oe Y) 0) 0 © do? 

oe Q, N oo F go ~ 

oe Su ~— Oo co O ce 

3° Gq Gg oe M% co Y 

oe “dam “AA oe oo TO 

oe ao — un — oO ft do © 
~— oe cd -d oe O ce eH 
qq oe ~~ ~~ oe oe QO 
© oe ame —~ bv de PH dO 
H oe “A “di oO oe 0) ce & 
Ex oe - Q ~~ N oe Mm dc? O 
oe dA vp~ 0 fx] do YD a 
= oe nn oO Nn co 12 ce O 
fd oe roKre) NO de fr] oe 4 
© oe — vo —- vO oe de 
fr, oe Y) n oe ch Yo 
OM oe om Oma oe co O 
4 oe Os oa oe de 

oe to er teed tees oe ce Oy 
we oe ~ ~ are oe oe O 
EX ae hoe wv oe oe oe Ha 

oe Doss DO es a_n oe 30° OY 
fr] oe ~~ Qy a ~ Ne eee fore) rave) C 
Ea oe Med aA eH~ ~~ ae oe YQ 
mG de eos G ~¢ dN oe de 
Zz oe ae eH OO ew HO & & axe) ce U 
HH oe ~ . oe OP -H oO 
A oe [hove § Ov | O&O tol oe oe 
| oe ~ ~ oe deo Y 
O o do N ~ GN oe oe ¢ 
© oe & & teas de de - 
UO oo M Oo oe de 

oe ie) Gdn oe for) 
de oe 4 08s ce de dO gO 


u4 


. 
/ 


‘ 


Bt ae 


sym(’u3(t)’) 


. 
, 


1) *m2 ( 
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. 
i 


Shh gad 
area 


. 
a 


rhoBuH_I + rhoPBd_I + rhoHTI 


‘ 


, 


rhoBuH*mi ( 
rhoPBd*mi ( 
diff(rho,t).’ 


f 


:k) 
{ul u2 Oj 
[el 0 0] 


{R 0 0] 
rhoPBd_I 
simplify ( 


[mb1 mb2 mb3 mb4,. mb5 mb6é mb7]; 


= mbt(1 
ul = sym(’ul(t)’); 
rhoBuH 

rhoBuH_I 

rho 

V(:,:,21) 


syms e1 R mbl mb2 mb3 mb4 mb5 mb6é mb7 


= sym(’u4(t)’); 


mbt 
mb 
rhoHI_I 
rhoPBd = 
end 





(UBk) 





2 
oO 
Potential energy or rotor blade 


fo) 
Ke) 
% 





(DBk) 


. 
‘ 


69 


UBk1 + UBk2; 
= 1/2*Czeta(i)*(diff(zet(i),t))*2 + 


[Vzetal Vzeta2 Vzeta3 Vzeta4 Vzeta5 Vzeta6 


[Czetal Czeta2 Czeta3 Czeta4 Czeta5 Czeta6 


Vzetat(1:k); 


[Kel Ke2 Ke3 Ke4 Ke5 Ke6 Ke7]; 
Czetat (1:k); 


Ke = Ket(1:k); 
[Kd1l Kd2 Kd3 Kd4 Kd5 Kd6 Kd7] 


Kd = Kdt(1:k); 
[Ks1l Ks2 Ks3 Ks4 Ks5 Ks6 Ks7] 
Ks = Kst(i:k); 


° 
, 


$UBk3 = 1/4*Ks(i)*signum(zet (i)-z) *(zet(i)%*24+z%*2- 


2*zet(i)*z)+1/4*Ks (i) *signum(zet (i)+z)*(-zet(1i)*2-z2*2- 


2° Zeer (il) ey rls 2 "Ke (i) eet (3) 2) ey 2 KS a) 22 


DBk (i) ( 
Vzeta(i)*(diff(zet(i),t))*2*abs(diff(zet(i),t)); 


end 


UBk (i) 
Dissapative function of rotor blades 


Syms Czetal Czeta2 Czeta3 Czeta4 Czeta5 Czeta6 Czeta7 
syms Vzetal Vzeta2 Vzeta3 Vzeta4 Vzeta5 Vzeta6 Vzeta7 





syms Kel Ke2 Ke3 Ke4 Ke5 Ke6 Ke7 
syms Kdl Kd2 Kd3 Kd4 Kd5 Kd6 Kd7 
syms Ksl Ks2 Ks3 Ks4 Ks5 Ks6 Ks7 z 
syms DB1 DB2 DB3 DB4 DB5 DB6 DB7 


Ket 
Kdt 
Kst 
for i 
end 
Czetat 
Czeta7]; 
Czeta 
Vzetat 
Vzeta7] 
Vzeta 


% 








oe de 
ae ce Pp 
oe 0° -r 
oe dep 
oe ce OG 
oe de O, 
oe ceo 6 
fo\ed de U} 
oe de WY) 
de OP -d 
oe oe A 
ae de 

ae de ™ 
food de 

foe) oe 
ae de® 1G 
oe de 1D 
de deo 
oo oe 

oe do DN 
oe de OD 
oe de 
oe oe OW 
oe de 
oe oe YO 
oe de 

co M oe 4H 
o& 1D & 
OO WT) dO ed 
oe oe 
Oo fy de 
wo O oe WU 
oe de pL 
ce HH co O 
de O de A 
o© A de 

oo fx] do ™ 
oe —& de 

o© fz} de —~ 
foe) de 
de de fH 
de de 

oe Ce 
de ae Om 
oe de 
oe de YO 
oe de 
oe ae 
oe oe 

oe co O 
oe oe -d 
oe de 4 
oo oe YO 
oe oe 
fale de eed 
oe dP RG 
ge ae 


dc od dO dP 


. 
‘ 


. 
, 


ff£(u2,t)%*2 
placements 


1 
1S 


dd 


1Ze 
MOTION USING LAGRANGE’S 
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. 
, 
Sard 
Fa 

x 


and zet(2) 


. 
f 


. 
a 
. 
af 


= u7; 


‘f 


or zet(l1) 
DOFB] 


F9 
diff (DOF, t) 


[Fl F2 u]; 


a 
. 


, 
° 
‘ 


(DH) 
diff (dDOF,t) 


{ul,u2] 
[DOFF, 


1/2*M (1) *di ft (ul,.b)*2 + 172*M(2)*d 


LOK CL) tale AY 2 PK (2). ue 


DES 2727601) Par et (ul, e)) 2+ 
ful u2 u3 u4 ud u6 u7] 


[cl ¢c2] 


[M1 M2] 
[Kl K2] 
F8 = ué; 


Generalized forces on general 


DERIVATION OF EQUATIONS 0O 


EQUATION 


1/2*6(2) * (AL fF (u2,t) )°2+L/24*v (1) * (aff (ul, &) )*2*abs (ditt (ul 


pt) )ai/2*v (2) (ar fi(u2, €))*2*abs (di fl (u2,t)) 


syms Kl K2 cl c2 vl v2 M1 M2 
syms gl q2 q3 q4 q5 q6 q7 q8 q9 


syms Fl F2 F3 F4 F5 u3 


function 
K 

TP 

UF 

z 

ut 

Fl 
u5; 

F 
DOFF 
DOF 
aDOF 
adDOF 





(op) 
ey 
se) 
ioe) 
ey 
8) 
~ 
ey 
fe) 
\O 
5! 
Ke) 
wn 
fey 
Ko) 
st 
ey 
0 
om 
5 
so) 
N 
ey 
Ko) 
ie 
oF 
so) 
; 
7] 


syms ddql ddq2 ddq3 ddq4 ddq5 ddqé ddq7 ddq8 ddq9 


f 


[ ql q2 q3 g4 q5 q6 q7 g8 QQ] 


DOFdt = 





f 


n) 
[ dql dq2 dq3 dq4 dq5 dq6 dq7 dq8 dq9] 


DOFg = DOFgt (1 


dDOFqt = 





f 





f 


n) 
[ A€dqi ddq2 ddq3 ddq4 ddq5 ddq6 ddq7 ddq8 ddq9J] 


dDOFq = dDOFqt (1 


ddDOFqt = 


. 
f 


Substution set construction 


% 


ee. 


set2 
for i 


size of DOFg vector 
[setl maple(‘‘=*‘’,DOF(i),DOFq(i) ) 


[setl maple(’*‘=*’,dDOF(i),dDOFq(i)) 
[setl maple(’‘='’,ddDOF(i),ddDOFq(i) ) 


[set2 maple(’' 


° 
© 


n 


=1 


ie 


setl 
setl 
setl 
set2 
set2 
set2 


. 
, 


] 


i? 


, 


] 


‘*, DOPG (1) «DOF (2).) 


7 


] 


[set2 maple(’‘=*',dDOFq(i) , DOF (i) ) 
[set2 maple(’‘=*'’,ddDOFq(i),ddaDOF (i) ) 


1 


end 


. 
, 


maple(’convert’,set1,'’set’) 


setl 
set2 


. 
f 


maple(’convert’,set2,’set’) 


for i 


Pauper. 


U = U+UBk (1); 


T = T+TBk ( 


D1+DBk (i) ; 


D1 


end 


Switch of time dependent terms 


2 
6 


maple(‘subs’,set1,T); 


Temp 


Complete EOM 


° 
6 
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for i= 1i:n 
temp1 = diff (Temp, GDOFq(i)); 
temp2 maple(’subs’,set2,temp1); 
temp3 = diff(temp2,t); 
L1 = maple(’subs’,setl1,temp3); 
L2 = diff(Temp,DOFq(i)); 
tempL3 = maple(’subs’,set1,U); 
L3 = diff(tempL3,DOFq(i)); %* check dimensions of DOFg 
tempL4 = maple(’subs’,set1,D1); 
L4 = diff(tempL4,dDOFq(i)); % check dimension of DOFq 
EOM(i) = simplify (L1-L2+L3+L4-F(i)); % check dimension 
of EOM and F 


end 
EEEEEEEEEEEEEESEEEETEEEEEEEEEEEEEEEEEEEEEEESESEEEEEESESESESES 
% Elimination of excess baggage 
CESEECESCEEEEEESESEEEEEEEESESEEEEEEEEEEEETEESESESEEEEESESESESS 
tempSmin = [ signuml(1,q3-z) signuml(1,q4-z) signuml1(1,q5- 
z) signuml(1,q6-z) signuml(1,q7-z) signuml1(1,q8-z) 


signumi(1,q9-z) ]; 

tempSpls = [ signuml1(1,q3+z) signuml(1,q4+z) 

signuml (1,q5+z) signuml(1,q6+z) signuml (1,q7+z) 
signuml(1,q8+z) signuml(1,q9+z) J; 

tempSabs = [ maple(’abs’,1,dql) maple(’abs’,1,dq2) 
maple(’abs’,1,dq3) maple(’abs’,1,dq4) maple(’abs’,1,dq5) 
maple(’abs’,1,dq6) maple(’abs’,1,dq7) maple('’abs’,1,dq8) 
maple(‘abs’,1,dq9)]; 

tempS = [tempSmin(1:n) tempSpls(i:n) tempSabs(1:n) J]; 
[g,p] = size(tempS) ; 

temp0O = zeros(1,p); 

setS = []; 


{setS maple(’‘=‘'‘,tempS(i),temp0(i)) ]; 


LESEESEEESEEEESEEESESESEEESEEESEEEEEEESEESESEEESEESEEEEEESESS 
% GENERATE Al AND £1 FROM EOM 
EESSESSEEEEEEEEEEESESEEEEEEEEEEEEEEEEESESESESEEEESEEEEESESS 
% Determine second derivative terms 
ESESETESESESESEEEEEEEEESEESEEEESEESEESEESESESEEESEESEEEEESESSEES 
for i = 1l:n 

for j = 1:n 

A(i,j) = maple(‘coeff’,EOM(i),ddDOFqg(j)); 
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° 
i 


13)) 


maple(’subs’,setS,A(i 


A(i,j) 


end 


end 


= ddDOFq*A; 


Ax2dot 


ate second der 


tive terms.from EOM 


iva 


oe 


. 
a 


(EOM (i) -Ax2dot (i) ) 
maple(’subs’,setsS, f ( 


)); 


i 


(2) 


Laity CEC) yyy 


imp 


-(s 


) 


i 


( 


£ 


end 


Generate state vector 


°. 
Ke) 


syms xl x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 


. 
a 


[xl x2 x3 x4 x5] 


X1ldot 


. 
, 


[x6 x7 x8 x9 x10] 


[Xltemp Xidot] 


Xltemp = 


. 
7 


X1 


[DOFq ADOFq] ; 


DOFqtemp = 


Convert terms to state vector form 


2. 
x) 


1=1:2*n 
setx 


for 


. 
, 


] 


[setX maple(’*‘=*‘’,DOFqtemp(i),X1(i)) 


end 


maple(‘convert’,setX,'set’) 


f 


setx 


a 


maple(’subs’,setX,A) 
fil = maple(’subs’,setxX,f£); 


Al 


Pages 


from input 


ee 


Eliminate variable 


2 
Kes 


syms ul u2 u3 


‘, 


[ul u2 u3] 
Al = subs(A1,u,ut) 


f1 = subs(f1,u,ut) 


ut = 


° 


? 


i i, 
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APPENDIX E. MATLAB® WORKSHEET FOR COMPLEX (STRAUB) MODEL 
WITH VARIABLE NUMBER OF ROTOR BLADES 


function [Al, £1] = compmod() 


clear 
cle 


k = 3; % number of blades, up to 7 
syms Phil Phi2 Phi3 Phi4 Phi5 Phi6é Phi7 Omega t alpha 


tempPhi = [Phil Phi2 Phi3 Phi4 Phi5 Phi6é Phi7]; 
Phi = tempPhi(1:k); 
psi = Omega*t + Phi; 


Yl oS sym rie)" )\Ss 22 = -syml'e2tey*)s 

r= iPr re] 

Tlr = [1 0 0; 0 cos(r1) sin(r1); 0 -sin(r1) cos(r1)]; 
T2r = [cos(r2) 0 sin(r2); 01 0; -sin(r2) 0 cos(r2)]; 
mil = (Tlr*T2r).’ 


for i = 1:k 


m2(:,:,1) = [cos(psi(i)) sin(psi(i)) 0; -sin(psi(i)) 
cos(psi(i)) 0; 0 0 ij; 
end 
zetl = sym(‘'zetl(t)‘'); zet2 = sym('zet2(t)'); zet3 = 
sym(‘’zet3(t)’); zet4 = sym(’zet4(t)’); 
zet5 = sym(’zet6(t)’); zet6 = sym(’zet6(t)’); zet7 = 


sym(’zet7(t)’); 
tempzet = [ zetl zet2 zet3 zet4 zet5 zet6 zet7 }; 
zet = tempzet(1:k); 


betl = sym(’bet1(t)’); bet2 = sym(‘bet2(t)’); bet3 = 
sym(‘bet3(t)’‘); bet4 = sym(’bet4(t)’); 
bet5 = sym(’bet6(t)’); bet6 = sym(’bet6(t)’); bet7 = 


sym(’bet7(t)’); 
tempbet = [ betl bet2 bet3 bet4 bet5 bet6 bet7 ]; 
bet = tempbet (1:k); 


for i= 1:k 
T3Z(:,:,i 


ty = [cos(zet(i)) sin(zet(i)) 0; -sin(zet(i)) 
cos (zet(i)) 


) 
QO; 00 1]; 
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T2b(:,:,1) = [cos(bet(i)) 0 sin(bet(i)); 01 0; - 
sin(bet(i)) 0 cos(bet(i))]; 

T3p(:,:,1) = [cos(psi(i)) sin(psi(i)) 0; -sin(psi(i)) 
cos(psi(i)) 0; 0 0 1]; 
end 


For. a. = (Lek 
m3(:,:,1) = (T3z(:,:,1)*T2b(:,:,1)).'; 
m4(:,:,1) = 
(T3z(:,:,1)*T2b(:,:,1)*T3p(:,:,1)*Tlr*T2r).’; 
end 


syms h el R mbl mb2 mb3 mb4 mb5 mb6 mb7 


tempmb = [mbl mb2 mb3 mb4 mb5 mb6 mb7]; 
mb = tempmb(1:k); 
1 sym? wl{t)); 
u2 sym(’u2(t)’); 


rhoFI_I = [ul u2 0]; 
rhoHF = [0 0h]; 

rhoBuH = [el 0 0]; % changed from vector form 
rhoPBd = [R 0 0]; 

rhoHF_I = (ml*rhoHF.’).'; 


for 1 = 1:k 
rhoBuH_I(:,:,1) 
rhoPBa_I(:,:,1) 
Phos, 354) = 

(rhoFI_I+rhoHF_I+rhoBuH_I(:,:,i1)+rhoPBG_I(:,:,1)).’; 
V(:,7:,1) = Giff (rho(:,:,1),t); 

Vsar (sy ty) = V1 ey hyper! oy V2 yn VSS 5 2) 39 
TBk(:,:,1) = 1/2*mb(i)*Vsqr(:,:,1); 


m1*m2(:,:,i)*rhoBuH. ’ 
ml*m2(:,:,1)*m3(:,:,1)*rhoPBd.’; 


end 
EESEESSEESESESSEEESSESEESEEEESEESSESEEEESSEEEEEEEEESEEEEESSEEEEE 
EESSESECESSCSCEESESESS% 

% Potential Energy of kth blade (UBk) 
BESEESEEEESESEEESESESSESESEEEESSESEESEEESEEEEEESESEESEEEESESEEESS 
EESEESEEECESEEEEEEESS 


syms Kel Ke2 Ke3 Kdl Kd2 Kd3 Ksl1 Ks2 Ks3 z 
syms kf11 kf12 kf£13 kf14 k£15 kf16 kf17 kf31 kf32 kf33 kf34 
k£35 k£36 k£37 
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syms k1l11 k112 k113 k114 k115 k116 k117 k131 k132 k133 k134 
k135 k136 k137 

syms Ke4 Ke5 Ke6 Ke7 Kd4 Kd5 Kd6 Kd7 Ks4 Ks5 Ks6 Ks7 
Ket = [Kel Ke2 Ke3 Ke4 Ke5 Ke6 Ke7]; 

Ke = Ket(1:k); 

Kdt = [Kdl Kd2 Kd3 Kd4 Kd5 Kd6 Kd7]; 

Kd = Kdt(1:k); 

Kst = [Ksl Ks2 Ks3 Ks4 Ks5 Ks6 Ks7]; 

Ks = Kst(1:k); 

kf£1t = [kf11 k£12 kf13 kf14 k£15 kf16 kf17]; 

k£1 = kf1t(1:k); 

kf£3t = [kf£31 k£32 k£33 kf34 kf£35 kf36 k£37); 

k£3 = k£3t(1:k); 

klit = [k111 k112 k113 k114 k115 k116 k117]; 

k11 = k11t(1:3); 

k13t = [k131 k132 k133 k134 k135 k136 k137]; 

k13 = k13t(1:3); 


for i=1:3 

UBk1 = 1/2* (bet (i)*2*kf£1(i)+zet(i)*2*k11(i)); 

UBk3 1/4* (bet (i) *4*k£3 (i) +zet(i)*4*k13(i)); 
UBk(i) = UBk1 + UBk3; % + UBk3; not required for simple 
model 


end 
ESSSESEESSESESSESESESEEESEETEEEEEESEESSESEEEEEEEEEEEEEEEEEEESESS 
ESESEESESEEESESSESS 

% Dissapative Energy 
SESEESESEEEEESEEEEESEEEEEEESEEEEESEEEESESEESEESEEEEEEEESSSSES 
SESEEEESEEEEESESESSES 


syms DB1 DB2 DB3 Czetal Czeta2 Czeta3 Vzetal Vzeta2 Vzeta3 
syms cfl cf2 cf£3 cf4 cf5 cf6 cf7 

syms cll cl2 c13 ¢cl4 ¢c15 cl6 cl7 

cft = [cf1l cf2 cf3 cf4 cf£5 cf£6 cf7]; 

ef = cft(1:k); 

elt = [cll cl2 cl13 cl14 ¢cl15 ¢l16 cl17]; 

el = clt(1:3); 

Czeta = [Czetal Czeta2 Czeta3]; % convert to 7 blades 
matrix 

Vzeta = [Vzetal Vzeta2 Vzeta3]; 


23 


= 1 
(1) AAA27 GLEE (bet (i), 6) °2* eb (i) + 
L/ 202 ff (vera) 7t)*2*el (3) + 
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££ (r2,t)*2*CR2 
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(DF) 
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(TP) 
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. 
, 
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1/2*ul*2*KT1 + 1/2*u2*°2*KT2 
1/2*r1*2*KR1 + 1/2*r2*2*KR2 
tion Energy of Fuselage 
1/2*diff(ul,t)*2*CT1 + 1/2*diff(u2,t)*2*CT2 
L724 dL thr 0)" 2* CR L/ 24a 


TFt + TFr 
UFt + UFr 


TR = 1/2* (Qi ff (ul, t)*24*M(1) +d ff (u2,t)*2*M(2) )+ 
Dissapa 


Tre = /2t (ar et (eit) Om 2e ld se areit (2, t) e241 22— 


syms Kl K2 cl c2 vl v2 M1 M2 I11 122 112 
2° GLEE (rb) Fai fb (22,2) 4212) 3 


syms CT1 CT2 CR1 CR2 VT1 VT2 VR1 VR2 


% Kinetic Energy of Fuselage 
% Potential Energy of Fuselage 


syms KT1 KT2 KR1 KR2 


TF 
UFt 
UFr 
UF 
DFtv 
DFrv 





DFth = 


T/2*di ff (ul, t)“2*abs (diff (ul, ¢)) *VP1l+1/2* diff (u2,t)*2*absita 


1/2*diff(r1,t)*2*abs (diff (r1,t))*VR1+1/2*diff(r2,t)*2*abs(d 


EEE(r2, tC) ) *VR2; 


DFtv + DFrv + DFth + DFrh; 


DF 


oe 
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(Generalized Aerodynamic Forces) 


% Aerodynamic 


oe 
oe 
oe 
0° 
oe 
oe 
oe 
oe 
oe 
oe 
de 
oe 
oP 
oe 
oe 
oe 
oe 
oe 
oe 
oe 
oe 


syms thetl thet2 thet3 thet4 thet5 thet6 thet7 cdO mu 


lambda rhol ac 


thett 


[thet1l thet2 thet3 thet4 thet5 thet6 thet7]; 


do for 7 blades 


° 
6 


Omega*lambda*R]; 


{mu*Omega*R, 0, 


Vair 
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2gby?e 
:,1)*°2+UT(: 
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t(2,: 


tf 


(3, 


= sgqrt(UP(: 


-V_Bd 
-V_Bd_t 


Aoced -d -d 


rag i 


tpi) /ULCs 


ane 


thet (i) -atan (UP ( 


2,1) 


1/2*rhol*a*c* (aoa(: 


f 


dFbeta(: 


:,1)- 


, 


2,1) *UT(: 


7 


:,i)*UU(: 


7 


‘ 


cd0/a*UU(: 


2,1) 


1/2*rhol*a*c* (aoa(i)*UU(: 


Fi 


dFzeta(: 


:,i1)*UT 


:,1)+cd0/a*UU(:, 


ba 


:,1)*UP(: 


i 


Mbeta_k 


$1) *Cos( set (1 ))> 


a, 


= 0.7*R*2*dFbeta(: 


Mest gad 


Mzeta(: 


end 


= 0.7*R*2*dFzeta(:,:,i); 
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%® Derivation of Equations of Motion by Lagrangian Method 
BEEESSEEEESEEEEEEEEESEEEEEEEEESEEEEEEEEESEESEEEEESEEEESEESEEES 
EESEBSEEEESEEEESESEES 

syms 


DOFF = [ul, u2, rl, r2]; 
DOFB = []; 
DOFB = [bet,zet]; 


DOF = [DOFF, DOFB]; 
GDOF = diff (DOF,t); 
GdDOF = diff(dDOF,t); 
{l,n] = size(DOF) ; 


EESEESEEEEEESEESESEEEEELESEEEBEEESEEEEESEBEEESEESEEEEESSESEESS 
EESEEEEEEEEEEEESESEEES 

% Transformation between time dependent and independent 
notation by substituting sets 
EESSESEESESEEEEEEEEEETEESESEESESEEESESSSSESESESSSSSSESSEEESSSS 
EESESESESESEEEEEEEESESS 

syms ql q2 q3 q4 q5 q6 q7 q8 q9 ql0 dql dq2 dq3 dq4 dgq5 dqé 


dq7 dq8 dq9 dq10 

syms ddqi ddq2 ddq3 ddq4 ddq5 ddq6é ddq7 ddq8 ddq9 ddq10 
syms ql1 ql2 ql3 ql4 qi5 ql6 qi7 ql18 dqil dql2 dqi3 dql14 
dqi5 dql6é dqi7 dqi8 

syms ddqil ddql2 ddql13 ddql4 ddq1l5 ddql6 ddql17 ddq18 


DOFqt = [ql q2 q3 q4 q5 g6 q7 g8 g9 ql0 qil qil2 ql3 ql4 qi5 
qié ql7 q18]; 

DOFgq = DOFqt(1:n); 

ADOFgt = [dql dq2 dq3 dq4 dq5 dq6é dq7 dq8 dq9 dqi0 dqill 
dq1i2 dgqi3 dqli4 dqi5 dql6 dqi7 dqi8]; 

ADOFq = dDOFqt (1:n); 

dadDOFgqt = [ddql ddq2 ddq3 ddq4 ddq5 ddqé ddq7 ddq8 ddq9 
ddqi0 ddqli ddqi2 ddql3 ddqi4 ddqi5 ddql6é ddql17 ddq18]; 
adDOFq = ddDOFgt (i:n); 


setl = []; 
set2 = []; 
for i=l:n % size of DOF vector 
setl = [setl maple(’‘=*'’,DOF(i),DOFq(i)) 7 


] 
setl = [setl maple(’‘='’,d@DOF(i),dDOFq(i)) ]; 
setl = [setl maple(’‘='’,ddDOF(i),ddDOFq(i)) ]; 
] 
) 


Z 


set2 = [set2 maple(’‘=*'’,DOFq(i),DOF(i)) 


set2 = [set2 maple(’‘=*'’,daDOFq(i),dDOF(i)) J; 
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set2 = [set2 maple(’‘=‘',ddDOFq(i),ddDOF(i)) ]; 
end 


setl = maple('‘convert’,setl,’set’); 
set2 = maple(’convert’,set2,’set’); 


T = TF; 
U = UF; 
D1 = DF; 


GF = [0,0,0,0]; 


[1,n] = size(DOFB) ; 
for 2. = sk 

T = T+TBk(:,:,1); 

U = U+UBk(i); 

D1 = D1+DBk(i); 

GF = [GF Mbeta_k(:,:,i) Mzeta(:,:,i)]; 
end 


Temp = maple(‘subs’,set1,T); %Temp = subs(T,dDOF,dDOFq) ; 
save for time comparison 

GFq = maple(’subs’,set1,GF); %GFq = subs(GF,dDOF,dDOFq) ; 
save for time comparison : 


fim] size(DOF) ; 
for i= i1i:n 
tempi = diff(Temp, dDOFq(i)); 
temp2 = maple(’subs’,set2,temp1); 
temp3 = diff(temp2,t) ; 
L1 = maple(’subs’,set1,temp3); 
L2 = diff (Temp, DOFq(i)); 
tempL3 = maple(’subs’,setl1,U); 
L3 = diff(tempL3,DOFq(i)); % check dimensions of DOFq 
tempL4 = maple(’subs’,set1,D1); 
L4 = diff(tempL4,dDOFq(i)); % check dimension of DOFq 


6 EOM(i) = simplify (L1-L2+L3+L4-GFq(i)); 

EOM(i) = (L1-L2+L3+L4-GFq(i)); % input to allow 
computation, gags on simplify(EOM) 
end 


% code will generate EOM, will not perform operations on 
EOM, ie. simplify, eval, etc. 
% remainder of code should be correct, could not evaluate 
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ore) 
oe 
ode 
er) 
oe 
oe 
oe 
oe 
oe 
oe 
oe 


Format equations of motion into form A d2x/dt2 = £ 


° 
ae) 


maple(‘coeff‘,EOM(i),ddDOFq(j)); 


Goll 
daa 
™ 
Hoos 
“d 
“Tl z 
My O 
ce) G 
a w 
ze) 
Gq 
v 


= size(ddDOFq) ; 
zeros(1,g); 


[5a] 
setz 


-~eval (subs (EOM(1:g),setZ,ddDOFq) ) ; 


El sg) 


= ddDOFg*A; 


Ax2dot 


for i 


. 
, 


£ (i) 


(EOM (i) -Ax2dot (i) ) 


end 


Change notation for input into S-function 


9 
6 


syms x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 


x17 x18 


[xl x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 


x15 x16 x17 x18]; 


Xltemp = 


2*n); 


[Xltemp Xidot] 


X1itemp(n+t+1: 


X1 


- 
a 


X1 


[DOFq ADOFq] ; 


DOFqtemp 


:2*n 


=1 


for i 


° 
ta 


] 


‘' ,DOFqtemp(i),X1(i)) 


[setX maple(‘‘= 


setx 


end 


maple(’convert’,setX,’set’); 


setx 


‘ 


abs(1,x2) abs(1,x3) abs(1,x4)] 


= [abs(1,x1) 


tempsetxX1l 
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tempzeros = zeros(size(tempsetxX1)); 


Al = maple(’subs’,setX,A); 
f1 = maple(’subs’,setX,£); 
£2 = subs(f1,tempsetX1,tempzeros); 


syms ul u2 u3 

ut = [ul u2 u3]; 

Al = subs(Al1,u,ut); 
£2 = subs(f2,u,ut); 
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"APPENDIX F. MATLAB® FUNCTION TO ACCESS THE MAPLE FUNCTION 


‘ABS’ TO CALCULATE THE ABSOLUTE VALUE 


oe 


abs.m 


de 


This function is designed to call the Maple function 


fol) 


b and c and return the result 


‘abs’ with two inpust, 


foes 


abs (b,c) 


function a 


a = maple(’abs’,b,c); 
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‘SIGNUM”’ 


This function is designed to call the Maple function 


APPENDIX G. MATLAB® FUNCTION TO ACCESS THE MAPLE FUNCTION 
signum.m 





a 


evaluating one expression, 


‘signum’ 


oe 


87 


t 


signum(a) 


maple(‘signum’,a) 





function s 


s 
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APPENDIX H. MATLAB® FUNCTION TO ACCESS THE MAPLE FUNCTION 
‘SIGNUM’, TAKING THE DERIVATIVE 


$signuml .m 


ction 


n 


% This function is design to call the Maple fu 


using a variation of the function with 


‘signum’ 


fo 


signum(a,b) 


function s 


maple(’signum’,a,b); 


s = 
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APPENDIX I. MATLAB® FUNCTION MUNION, TAKING THE UNION OF 
TWO SETS USING THE MAPLE FUNCTION 


function sym_union = munion(argl1,arg2) 


stringa = char(argl); *%get rid of extra char’s 
stringa = stringa(8:end-2); 


if length(stringa) == 0 
stringa = ‘{}’; 
end 


stringb char (arg2); 
stringb = stringb(8:end-2); 


if length(stringb) == 0 

stringb = ‘{}’'; 
end 
union_string = [stringa,’ union ',stringb]; 
union_string (unionstring == *%]*) = °} "ys 
Wnion. String (union.string ==-" [ep =: "hss 
union_string = [’maple(’’ ‘, union_string, rE yey s 
temp1 = sym(eval (union_string) ); 
temp2 = char(templ); %@[templ(templ == ’{’) = ’[’] 


ta _ Ca 


temp2 (temp2° ==" {*) = 71 
temp2(temp2 == '}’) = '] 


ta 


, 
foe 
vy 


sym_union = sym(temp2) ; 
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